{% load static %}
library(sf) # simple features for geometries
library(rgdal) # OGR GeoJSON driver for importing data
library(leaflet) # mapping visualization
url <- "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month.geojson" # read geojson url into a variable
earthquakes <- readOGR(url, verbose=FALSE) # use the rgdal library to read the url into a SpatialPointsDataFrame
eqsf <- st_as_sf(earthquakes) # convert sp data into a simple features dataframe
# use the query tool in the REST page (set where clause to '1=1' to get all results) to generate a geojson request url
faults_url <- "https://services.arcgis.com/jIL9msH9OI208GCb/ArcGIS/rest/services/Active_Faults/FeatureServer/0/query?where=1%3D1&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&resultType=none&distance=0.0&units=esriSRUnit_Meter&relationParam=&returnGeodetic=false&outFields=*&returnGeometry=true&featureEncoding=esriDefault&multipatchOption=xyFootprint&maxAllowableOffset=&geometryPrecision=&outSR=&defaultSR=&datumTransformation=&applyVCSProjection=false&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnExtentOnly=false&returnQueryGeometry=false&returnDistinctValues=false&cacheHint=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&returnZ=false&returnM=false&returnExceededLimitFeatures=true&quantizationParameters=&sqlFormat=none&f=pgeojson&token="
faults <-readOGR(faults_url, verbose=FALSE) # use the rgdal library to read the url into a SpatialPointsDataFrame
faults_sf <- st_as_sf(faults) # convert sp data into a simple features dataframe
# create class breaks with "Spectral" based on earthquake magnitude
pal <- colorBin(
palette = "Spectral",
domain = eqsf$mag, # use magnitude variable
reverse = TRUE, # reverse the color direction
bins = 5 # 5 breaks
)
leaflet() %>%
setView(-117.841293, 46.195042, 3) %>% # set view to greater North America
addProviderTiles(providers$CartoDB, group = "Grayscale") %>% # add CARTO tiles
addProviderTiles(providers$Esri.WorldTerrain, group = "Terrain") %>% # add esri tiles
addPolylines(data = faults_sf,
popup = paste0( # create custom pop up
"<strong>Name:</strong> ", faults_sf$name,
"<br><strong>Slip Type:</strong> ", faults_sf$slip_type
),
group = "vectorData") %>% # add fault lines
addCircleMarkers(data = eqsf, # create circle markers from earthquake data
fillColor = ~pal(mag),
radius = ~eqsf$mag * 2,
stroke = FALSE,
color = "Spectral",
fillOpacity = 0.6,
popup = paste0( # create custum pop-ups
"<strong>Title:</strong> ", eqsf$title,
"<br><strong>Magnitude:</strong> ", eqsf$mag,
"<br><strong>Intensity:</strong> ", eqsf$mmi,
"<br><strong>Significanceq:</strong> ", eqsf$sig
),
group = "vectorData") %>%
# add legend to the bottom left
addLegend(pal = pal, values = eqsf$mag, position = "bottomleft", title = "Magnitude") %>%
# create a layer toggle for the basemaps and vector data
addLayersControl(overlayGroups = c("vectorData"), baseGroups = c("Terrain", "Grayscale"))